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Abstract. The present study derives the distribution of secondary masses M2 for the 67 exoplanets and very 
low-mass brown-dwarf companions of solar-type stars, known as of April 4, 2001. This distribution is related to 
the distribution of M2 sin i through an integral equation of Abel's type. Although a formal solution exists for 
this equation, it is known to be ill-behaved, and thus very sensitive to the statistical noise present in the input 
M2 sin i distribution. To overcome that difficulty, we present two robust, independent approaches: (i) the formal 
solution of the integral equation is numerically computed after performing an optimal smoothing of the input 
distribution, (ii) the Lucy-Richardson algorithm is used to invert the integral equation. Both approaches give 
consistent results. The resulting statistical distribution of exoplanet true masses reveals that there is no reason 
to ascribe the transition between giant planets and brown dwarfs to the threshold mass for Deuterium ignition 
(about 13.6 Mj). The M2 distribution shows instead that most of the objects have M2 < 10 Mj, but there is a 
small tail with a few heavier candidates around 15 Mj. 
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1. Introduction 



Han et al. (2001) suggested that most of the exoplanet 
candidates discovered so far have masses well above the 
lower limit defined by sinz = 1 (where i is the inclina- 
tion of the orbital plane on the sky) and should therefore 
be considered as brown dwarfs or even stars rather than 



planets . The present paper shows that this claim is not 
consistent with the distribution of masses extracted from 
the observed M2 sin i distribution (where M% is the mass 
of the companion) under the reasonable aooumption that 



The basic reason why the M% distribution obtained 
in Sect. [| differs from that of Han et ah (2001) is be- 
cause these authors concluded that most of the sys- 
tems containing exoplanet candidates are seen nearly 
pole-on. This conclusion, based on the analysis of the 
Hipparcos Intermediate Astrometric Data (IAD), has how- 
ever been shown to be incorrect (e.g. Halbwachs et al 



2000| ; |Pourbaix 200f ; |Pourbaix fc Arenou 2001[) , as sum- 
marized in Sect. ^[ 



While this paper was being referred, 



Zucker & Mazeh 



the ort Mto are oriented at random in opaco. Although the 
distributions of Mi and M2 sin i are related through an in- 



(2001| ) and [Tabachnik fc Tr emainc (200l|) proposed other 



tegral equation of Abel's kind (Chandrasekhar fc Munch 



195C; Lucy 1974), its numerical solution is ill-behaved. 



Two different methods are used here to overcome that dif- 
ficulty. In the first method (Sect. ||), the formal solution of 
Abel's equation is implemented numerically on an input 
M2 sin i distribution that has been optimally smoothed 



with an adaptive kernel procedure (Silverman 1986) to 
remove the high-frequency fluctuations caused by small- 
number statistics. The other method (Sect. p|) is based 
on the Lucy-Richardson inversion algorithm ( [Richardson 
1972|; [Lucy 1974]). 
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interesting approaches to derive the exoplanet mass dis- 
tribution. Zucker & Mazeh derive the binned true mass 
distribution by using a maximum likelihood approach to 
retrieve the average values of the mass distribution over 
the selected bins. Their results are in very good agree- 
ment with ours. Tabachnik & Tremaine suppose that the 
period and mass distributions follow power laws, and de- 
rive the corresponding power-law indices from a maximum 
likelihood method. On the contrary, the methods used in 
the present paper (and in Zucker & Mazeh's) are non- 
parametric in nature, since they do not require to make 
any a priori assumptions on the functional form of the 
mass distribution. This is especially important since the 
comparison of the shapes of the mass distributions for 
exoplanets and low-mass stellar companions may provide 
clues to identify the process by which they formed. By 
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imposing a power-law function like that observed for low- 
mass stellar companions, Tabachnik & Trcmainc (2001) 
somehow implicitly assume that these processes must be 
similar. 



2. The integral equation of Abel's kind relating 
the distributions of M2 sini and Mi 

The M 2 sin i values for low- mass companions of main se- 
quence stars may be extracted from the spectroscopic 
mass function and from the primary mass as derived 
through e.g., isochrone fitting. Let $>(Y) be the observed 
distribution of Y = M2 sin i which is easily derived from 
the observed spectroscopic mass functions provided that 
M2 << Mi as it is expected to be the case for the systems 
under consideration. Then, the seeked distribution \t(M 2 ) 
obeys the relation 



*(M 2 ) T1(Y\M 2 ) dM 2 . 



(1) 



The kernel II(F|M 2 ) corresponds to the conditional prob- 
ability of observing the value Y given M 2 . Under the as- 
sumption that the orbits are oriented at random in space, 
the inclination angle i distributes as sini, and the follow- 
ing expression is obtained for the kernel: 



n(Y\M 2 ) 



sinjo 
M 2 cos io 1 



(2) 



where io satisfies the relations M 2 sin io — Y — and 
< io < 90. Eliminating the inclination io i n the above 
relation yields 



n(F|M 2 



Y 



1 



m 2 (M| - y 2 )V2 

Eq.(|l|) thus rewrites 

/oo 
*(M 2 



with Y < M 2 . (3) 



M 2 (M 2 -Y 2 y/ 2 



dM 9 



(4) 



Eq.(Q) is the integral equation to be solved for <3>(M 2 ). It 
can be reduced to Abel's integral equation by the substi- 
tutions ( phandrasckhar fc Munch 1950] ) 

Y 2 = I/77 and A/f = 1/f. (5) 

With these substitutions, Eq.ffl) becomes 



(v-t) 1 / 2 



where 



Hv) = *(-J=) and m = _L-tf (-L). 



(X) 



It is well known that the solution of Abel's equation 
(Eq. ||) is given by 



m = - 



Wo dr,{t~ V y/ 2 dV+ it y/i 



where </>(0) = lim Y _ >00 $(y) = 0. 

While Eq. (^) represents the formal solution of the 
problem, it is difficult to implement numerically, since it 
requires the differentiation of the observed frequency dis- 
tribution $>(Y). Unless the observations are of high pre- 
cision, it is well known that this process can lead to mis- 
leading results. To overcome that difficulty, the observed 
frequency distribution has been smoothed in an optimal 
way (see Appendix) before being used in Eq. (||). The so- 
lution ^(t) is then computed numerically using standard 
differentiation and integration schemes. 

3. The Lucy-Richardson inversion algorithm 
applied to Abel's integral equation 

The Lucy-Richardson algorithm provides another robust 
way to invert Eq. (g) (see also |Cerf fc Boffin 1994| ). The 
method starts from the Bayes theorem on conditional 
probability in the form 

tf(M a ) n(Y|M a ) = $fY) R(M 2 \Y), (9) 

where R{M 2 \Y) is the reciprocal kernel corresponding to 
the integral equation inverse to the one that needs to be 
solved (Eq. 0): 

*(M 2 )=/ R(M 2 \Y) dY, (10) 



The reciprocal kernel represents the conditional prob- 
ability that the binary system has a companion mass M 2 
when the observed M 2 sini value amounts to Y. Thus, one 
has: 

*(m 2 ) n(y|M 2 ) 



R(M 2 \Y) = 



*(M 2 ) U(Y\M 2 ) 



J °° *(Af 2 ) U{Y\M 2 ) dM 2 



(11) 
(12) 



which obviously satisfies the normalization condition 
J °° R{M 2 \Y) dM 2 = 1. The problem in solving Eq. (jToj) 
arises because R(M 2 \Y) also depends on <3>(M 2 ), so that 
an iterative procedure must be used. If ^> r (M 2 ) is the rth 
estimate of ^(M 2 ), it can be used to obtain the (r + l)th 
estimate in the following way: 

/■Ma 

tt r+ i(M 2 ) = / $(F) R r (M 2 \Y) dY (13) 



with 



(6) R r (M 2 \Y) = 



and 
$ r (F) 



VAM 2 ) U(Y\M 2 ) 



* r (M 2 ) n(F|M 2 ) dM 2 . 



(14) 



(15) 



Thus, & r (Y) represents the corresponding rth estimate 
of the observed distribution $(Y). Eqs.(|l3|) and (14) to- 
gether yield the recurrence relation for the 'Iv's, 

,M 2 $(y) 



(8) * r+1 (M 2 ) = tt r (M 2 ) 



MY) 



U(Y\M 2 ) dY, (16) 
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with n(Y|M 2 ) given by Eq. (||) for the problem under 
consideration. The conditions for convergence of this re- 
currence relation are discussed by [Lucy (1974 ) and 3erf & 



Boffin [1994 ). It needs only be remarked here that (i) the 
iterative scheme converges if & r (Y) tends to &(Y), given 
the normalization of the probability n(Y|M 2 )dy, and (ii) 
the full convergence of the method is not necessarily desir- 
able, as the successive estimates $ r (Y) will tend to match 
$>(Y) on increasingly smaller scales, but the small-scale 
structure in &(Y) is likely to be dominated by the noise 
in the input data. This is well illustrated in Fig. |2] below. 

When the number of data points is small (typically 
N < 100; Ccrf k Boffin 1994 ), it is advantageous to ex- 
press $(y) as 



1 N 



(17) 



where the y n (n = 1, . . . N) are the N individual measured 
Mi sini values and S(x) is the Dirac 'function' such that 
xq = J 5(x — xq) dx. Substitution in Eq. ( |T3| ) then yields 



1 N 

* r+ i(M 2 ) = -^i? r (M 2 



(18) 



where i? r (M 2 |y„) is defined as in Eq. (|14|). The sample 
size should nevertheless be large enough for the functions 
Rr(M<2\y n ) to have sufficient overlaps so as to produce a 
smooth ^ r+ i(M 2 ) function. 

In the application of the method described in Sect. |], 
the initial mass distribution tE^C^a) was taken as a uni- 
form distribution, but it has been checked that this choice 
has no influence on the final solution. 



4. The frequency distribution ^(il^) 

The cumulative frequency distribution of the M 2 sin i val- 
ues smaller than 17 Mj, where Mj is the mass of Jupiter 
(= M Q /1047.35), available in the literature (as of April 
4, 2001) is presented in Fig. [3. It appears to be sufficiently 
well sampled to attempt the inversion procedure. The 
corresponding frequency distributions smoothed with two 
different smoothing lengths, locally self-adapting around 
h pt = 1 Mj and 2h op t (see Appendix) are presented as 
well for comparison. 

The sample includes 60 main-sequence stars hosting 
67 companions with Af 2 sini < 17 Mj. Among those, 6 
stars are orbited by more than one companion, namely 
HD 74156 and HD 82943 (M 2 smi = 1.5 0, 7.40 Mj, and 
M 2 si n^ = 0.84, 1.57 Mj, respecti vely; |Udry fc Mayor 
200l])F|, HD 83433 (0.16, 0.38 Mj; |Mayor et al. 2000|) 



HD 168443 (7.22, 16.2Mj;[Udry et al. 2000b|) , v And (0.71, 
2.20 and 4.45 M, T ; [Butler et al. 1999[ ) and Gliese 876 (0.56, 
1.88 Mj; Marcy et al. 2001). The inversion process is only 



able to treat these systems under the hypothesis that the 



see also http://www.eso.org/outreach/press-rel/pr- 



200 l/pr-07-01 .html 
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Fig. 1. Bottom panel: The cumulative frequency distribu- 
tion of the observed M 2 sin i values for the 67 exoplanets 
and very low-mass brown dwarfs (around 60 stars) known 
as of April 4, 2001. Top panel: The corresponding fre- 
quency distribution, represented either as an histogram, 
or smoothed using an Epanechnikov kernel (see Appendix) 
with an optimal smoothing length of h opt = 1 Mj (dashed 
curve). The frequency distribution smoothed with 2/i opt is 
also shown for comparison (solid curve) 



orbits of the different planets in a given system are not 
coplanar, since Eq. (^) to hold requires random orbital in- 
clinations. The case of coplanar and non-coplanar orbits 
are discussed separately in the remaining of this section. 
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4.1. Non-coplanar orbits 

Figure || compares the solutions ^(il^) obtained from the 
Lucy-Richardson algorithm (after 2 and 20 iterations, de- 
noted *2 and ^20, respectively) and from the formal so- 
lution of Abel's integral equation with smoothing lengths 
h op t and 2 /i opt on &(Y) (the corresponding solutions are 
denoted ^>h opt and ^2h opt )- The solutions from the two 
methods basically agree with each other, although solu- 
tions with different degrees of smoothness may be ob- 
tained with each method. On the one hand, Vt^o and *h opt 
exhibit high-frequency fluctuations that may be traced 
back to the statistical fluctuations in the input data. This 
can be seen by noting that the peaks present in ^>h opt cor- 
respond in fact to the high-frequency fluctuations already 
present in <£fo opt (Fig. |l|). These fluctuations should thus 
not be given much credit. The same explanation holds true 
for <3>20, since it was argued in Sect. || that the solutions 
\l/ r . resulting from a large number of iterations tend to 
match $ at increasingly small scales (i.e., higher frequen- 
cies) where statistical fluctuations become dominant. On 
the other hand, ^2 and ^2h opt are much smoother, and 
are probably better matches to the actual distribution. 
The local maximum around M2 ~ 1 Mj is very likely 
however an artifact of the strong detection bias against 
low-mass companions. 

The most striking feature of the 4"(M 2 ) distribution 
displayed in Fig. || is its decreasing character, reaching 
zero for the first time around M2 = 10 Mj, and in any 
case well before 13.6 Mj. The latter value, correspond- 
ing to the minimum stellar mass for igniting Deuterium, 
does not in any way mark the transition between giant 
planets and brown dwarfs, as sometimes proposed. That 
transition, which is thus likely to occur at smaller masses, 
must rely instead on the different mechanisms governing 
the formation of planets and brown dwarfs. Another ar- 
gument favouring a giant-planet/brown-dwarf transition 
mass smaller than 13.6 Mj is provided e.g., by the ob- 
servation of free-floating (and thus most likely stellar) ob- 
jects with masses pro bably smaller than 10 Mj in the a 
Orionis star cluster ( [Zapatcro Osorio et al. 2000| ). The 
\&(M2) distribution nevertheless clearly exhibits a tail of 
objects clustering around M2 ~ 15 Mj, due to HD 114762 
(M 2 sini = 11.5 Mj), HD 162020 (14.3 Mj), HD 202206 
(15.0 Mj) and HD 168443c (16.2 Mj). It would be interest- 
ing to investigate whether these systems differ from those 
with smaller masses in some identifiable way (periods, ec- 
centricities, metallicities,...), so as to assess whether or not 
they form a distinct class (Udry et al., in preparation). 
The jackknife method (e.g., Lupton 1993] ) has been 



used to estimate the uncertainty on the ^2h opt solution. 
In a first step, 67 input &2h opt distributions are computed, 
corresponding to all 67 possible sets with one data point 
removed from the original set. Eq. || is then applied to 
these 67 different input distributions. The resulting dis- 
tributions are displayed in Fig. |[ which shows that the 
threshold observed at 10 Mj is a robust result not affected 
by the uncertainty on the solution. 



0.4 



0.3 - 



0.2 - 



0.1 




Fig. 2. Comparison of ^(A^) distributions obtained from 
the Lucy-Richardson algorithm (after 2 and 20 iterations, 
represented by the solid and dashed histograms, respec- 
tively), and from the formal solution of Abel's integral 
equation with smoothing lengths 2h opt and h opt applied 
on $(F) (represented by solid and dashed curves, respec- 
tively). The supposedly best representation of the actual 
^{M-i) distribution is represented by the thick solid line. 
Planetary systems have been included in the inversion pro- 
cess with the assumption of non-coplanarity of the plane- 
tary orbits in a given system 



4.2. Coplanar orbits in multi-planets systems 

All the results discussed so far are obtained under the as- 
sumption that orbits of planets belonging to a planetary 
system are not coplanar. To evaluate the impact of this 
hypothesis, the following procedure has been applied. In 
a first step, the Lucy-Richardson algorithm is applied on 
the data set excluding the 13 planets belonging to plan- 
etary systems. That mass distribution obtained after 2 
iterations is then completed by mass estimates for the re- 
maining 13 planets. For each of the 6 different systems, 
an inclination i is drawn from a sini distribution. This is 
done through the expression i — acos x, where a; is a ran- 
dom number with uniform deviate. The same value of i is 
then applied to all planets in a given system to extract M2 
from the observed M2 sini value. The distributions of ex- 
oplanet masses obtained with and without the hypothesis 
of coplanarity are compared in Fig. ^, and it is seen that 
planetary systems are not yet numerous enough for the 
coplanarity hypothesis to alter significantly the resulting 
*(M 2 ) distribution. 

In any case, the main result of the present paper is 
that the statistical properties of the observed M% sin i dis- 
tribution coupled with the hypothesis of randomly oriented 
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Fig. 3. Comparison of the input $2fi opt distribution (thick 
dashed line) and the 67 ^2h op t distributions (thin dot- 
ted lines) resulting from the application of the jackknife 
method (see text), which illustrates the uncertainty on 
the solution x &2h apt (thick solid line). To guide the eye, a 
power-law of index —1.6 has been plotted as well (dot- 
dashed line) 
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Fig. 4. Evaluation of the impact of the coplanarity hy- 
pothesis on the resulting ^(A'h) distribution. The dashed 
histogram corresponds to the mass distribution obtained 
assuming coplanar orbits in planetary systems (see text), 
as compared to the ^2 (solid histogram) and ^2h apt solu- 
tions (see Fig. ||) 



orbital planes confine the vast majority of planetary com- 
panion masses below about 10 Mj. Zuckcr & Mazeh (2001) 
reach the same conclusion. 



It should be remarked that the above conclusion can- 
not be due to detection biases, since the high-mass tail 
of the M2 distribution is not affected by the difficulty of 
finding low-amplitude, long-period orbits. 



5. Discussion 

Although the assumption of random orbital inclinations 
seems reasonable, it is at variance with the conclusion of 



Han ct al. (2001 ) that most of the systems containing ex- 
oplanet candidates are seen nearly pole-on. These authors 
reached this conclusion by trying to extract the astromct- 
ric orbit, hence the orbital inclination, from the Hipparcos 
IAD. Halbwachs et al. (2000 ) had already cautioned that 
this approach is doomed to fail for systems with apparent 
separations on the sky that are below the Hipparcos sensi- 
tivity (i.e. ^ 1 mas). In those cases, the solution retrieved 
from the fit of the IAD residuals is spurious, since the true 
angular semi-major axis a is simply too small to be seen by 
Hipparcos. Since Halbwachs ct al. (2000] ) have shown that 
a actually follows a Rayleigh probability distribution, the 
fit of the IAD residuals will yield a solution larger than the 
true value, in fact of the order of the residuals. But since 
a sin i is constrained by the spectroscopic orbital elements, 
the too-large astromctric a value will force i to be close to 
to match the spectroscopic value of the product asini, 



as convincingly shown by Pourbaix (2001). Hence, this ap- 
proach gives the impression that all orbits are seen nearly 
face-o n. As an illustrative example, |Pourbaix fc Arcnor 



(2001 ) have shown that such an approach leads to a stel- 



lar mass for the companion of HD 209458 that, on another 
hand, has been proven to be a 0.69 Mj planet by the pho- 
tometric observation of the planet transit in front of the 
star ( |Charbonncau ct al. 200C ). 



The Han et al. (2001 ) result is moreover statistically 



very unlikely if the orbital planes are oriented at random 
in space (Pourbaix & Arenou 2001). Han ct al. (200l| ) 



have tried to justify this unlikely statistical occurrence by 
invoking biases against high-amplitude orbits in the se- 
lection process of the radial-velocity-monitoring samples. 
To the contrary, the planet-search surveys were specifi- 
cally devised to avoid such biases, as they aim at find- 
ing not only giant planets but also brown dwarfs so as to 
constrain the substellar secondary mass function of solar- 
type stars. Furthermore, Han et al. (2001) argument is 



totally invalid in the case of volume-limited, statistically 
well-defined samples like that of the CORALIE planet- 



search programme in the southern hemisphere ( Udry et al 
|2000a|) . This sample has been specifically designed to 



detect companions of solar-type stars all the way from 
q = M2/A/1 = 1 down to q < 0.001. 
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Appendix: Non-parametric treatment of the data 

To decrease the noise and allow a tractable use of the in- 
formation present in small data samples, heavy smoothing 
techniques are often required. A common practice con- 
sists in converting a set of discrete positions into binned 
"counts" . Binning is a crude sort of smoothing and many 
studies in statistical analysis have shown that, unless the 
smoothing is done in an optimum way, some, or even most, 
of the information content of the data could be lost. This 
is especially true when a large amount of smoothing is 
necessary, which then changes the "shape" of the result- 
ing function. In statistical terms, the smoothing process 
not only decreases the noise (i.e., the function's variance), 
but at the same time introduces a bias in the estimate of 
the function. 

The variance-bias trade off is unavoidable but, for a 
given level of variance, the bias increase can be minimized. 
The correct manner of achieving that task is provided by 
the so-called non-parametric density estimate methods for 
the determination of the "frequency" function of a given 
parameter or by the non-parametric regression methods 
for the determination of a smooth function g inferred from 
direct measurements of g itself. Moreover, adaptive non- 
parametric methods are designed to filter the data in some 
local, objective way minimizing the bias, in order to get 
the smooth desired function without constraining its form 
in any way. The theory and algorithms related to those 
methods, originally built to handle ill-conditioned prob- 
lems (either under-determined or extremely sensitive to 
errors or incompleteness in the data), are widely discussed 
in the specialized literature and summarized in easy-to- 
read textbooks (e.g., Silverman 1986| ; Hardle 1990; Scott 
1992j)7| 

The simplest of the available algorithms is provided by 
the kernel estimator leading to the following form of the 
normalized "frequency" function 



AT 



71=1 V 

with the normalization 



K(u) du = 1, 



(19) 



(20) 



where X n (n = l,...N) are the n available data points. 
Each data point is thus simply replaced by a 'bump'. 
The kernel function K determines the shape of the bumps 
while the bandwidth h (also called smoothing parameter) 
determines their widtb^. 

2 In the regres s ion problem, the Nadaraya- Watson estimator 
( Nadaraya 1964| ; Watson 1964) is commonly used: 



E 

n.=l 



9(X n )K{- 



X 



-) 



In the adaptive kernel version, a local bandwidth h n = 
h{X n , /) is defined and used in Eq. [n| In order to follow 
the "true" underlying function in the best possible way, 
the amount of smoothing should be small when / is large 
whereas more smoothing is needed where / takes lower 
values. A convenient method to do so consists in deriving 
first a pilot estimate / of /, e.g. by using an histogram or 
a kernel with fixed bandwidth /i op t, and then by defining 
the local bandwidths 



h n = h(X n ) = h opt [f(X n )/s] 
where 

1 N 

lo gS ,-^log/(I„). 



(21) 



(22) 



n=l 



It may be shown (Abramson 1982) that a = 1/2 leads to 
a bias of smaller order than for the fixed h estimate. 

The optimum kernel K may be taken as the one min- 
imizing the integrated mean square error beween / and / 
(MISE), where 



MISE(f) = E / f(x) - f(x) 



dx 



(23) 



is usually taken as an indicator of efficiency. In the above 
expression, E denotes the statistical expectation. This 
optimum kernel is the so-called Epanechnikov (K e ) or 
quadratic kernel: 



K e {u) = -(l-u 2 ), \u\<l. 



(24) 



Other choices of K differ only slightly in their asymptotic 
efficiencies and could be more adapted to particular pur- 
poses as e.g., for bi-dimensional data. 

The pilot smoothing length (h op t) is the only subjec- 
tive parameter required by the method. It relates to the 
quality of the sampling of the variable under considera- 
tion. There are several ways for automatically estimating 
an optimum value of h opt (see e.g. Silverman 1986, for an 
extensive review). A simple approach based on the data 
variance gives in our case h op t — 1.0 Mj. As the derivative 
of the frequency function rather than the function itself is 
actually used in Eq. (g), a larger pilot smoothing length 
(2h op t) was also considered in order to remove spurious 
small-statistics fluctuations of the density estimate. 

References 

Abramson I. S. 1982, Ann. Statist., 10, 1217 
Butler R. P., Marcy G. W., Fischer D. A. et al. 1999, ApJ, 
526, 916 

Cerf N., Boffin H. M. J. 1994, Inverse Problems, 10, 533 
Chandrasekhar S., Munch G. 1950, ApJ, 111, 142 
Charbonneau D., Brown T. M., Latham D. W., Mayor M. 

2000, ApJ, 529, L45 
Halbwachs J. L., Arenou F., Mayor M., Udry S., Queloz 

D. 2000, A&A, 355, 581 
Han I., Black D., Gatewood G. 2001, ApJ, 548, L57 



A. Jorissen et al.: The distribution of exoplanet masses 



7 



Hardlc W. 1990, Applied nonparamctric regression, 

Cambridge University Press, Cambridge 
Lucy L. B. 1974, AJ, 79, 745 

Lupton R. 1993, Statistics in theory and practice, 

Princeton University Press, Princeton 
Marcy G. W., Butler R. P., Fischer D. A., Vogt S., Lissauer 

J. 2001, ApJ, 556, 296 
Mayor M., Nacf D., Pepc F. et al. 2000, in A. Penny, 
P. Artymowicz, A.-M. Lagrange, S. Russell (eds.), 
Planetary Systems in the Universe: Observations, 
Formation and Evolution, ASP Conf. Ser., in press 
Nadaraya E. A. 1964, Theor. Probab. Appl., 10, 186 
Pourbaix D. 2001, A&A, 369, L22 
Pourbaix D., Arcnou F. 2001, A&A, 372, 935 
Richardson W. H. 1972, J. Opt. Soc. Am., 62, 55 
Scott D. W. 1992, Multivariate Density Estimation, 

Multivariate Density Estimation, Wiley, New York 
Silverman B. W. 1986, Density estimation for statis- 
tics and data analysis, Monographs on Statistics and 
Applied Probability, London: Chapman and Hall 
Tabachnik S., Tremaine S. 2001, AJ, in press 
Udry S., Mayor M. 2001, in Astrobiology, Lecture Notes 

in Physics, Springer Verlag 
Udry S., Mayor M., Naef D., Pope F., Queloz D., Santos 
N. C, Burnet M., Confino B., Melo C. 2000a, A&A, 
356, 590 

Udry S., Mayor M., Queloz D. 2000b, in A. Penny, 
P. Artymowicz, A.-M. Lagrange, S. Russell (eds.), 
Planetary Systems in the Universe: Observations, 
Formation and Evolution, ASP Conf. Ser., in press 

Watson G. S. 1964, Sankhya, Series A, 26, 359 

Zapatero Osorio M. R., Bejar V. J. S., Martin E. L., 
Rcbolo R., Barrado Y Navascues D., Bailer- Jones C. 
A. L., Mundt R. 2000, Science, 290, 103 

Zucker S., Mazch T. 2001, ApJ, in press 



Mass distribution of exoplanets 
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